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Abstract. Dynamics of linear perturbations in a differen¬ 
tially rotating accretion disk with non-homogeneous ver¬ 
tical structure is investigated. It has been found that tur¬ 
bulent viscosity results in instability of both pinching os¬ 
cillations, and bending modes. Not only the low-frequency 
fundamental modes, but also the high-frequency reflective 
harmonics appear to be unstable. The question of the lim¬ 
its of applicability of the thin disk model (MTD) is also in¬ 
vestigated. The insignificant distinctions in the dispersion 
properties of MTD and three-dimensional model appear 
for wave numbers k < (1 -y 3 )/h (h is the half-thickness 
of a disk). In the long-wavelength limit, the relative dif¬ 
ference between eigenfrequencies of the unstable acoustic 
mode in the 3D-model and the MTD is smaller than 5 %. 

Key words: accretion disk — instability — turbulent 
viscosity — acoustic mode 


1. Introduction 

Accretion disks (AD) are important in many observable 
astrophysical objects (close binary systems, quasars, ac¬ 
tive galactic nuclei, young stars, protoplanet disks). It is 
usually assumed that the disk is geometrically thin, and 
that the disk has turbulent viscosity (Shakura & Sunyaev 
1973; Lightman & Eardley 1974). The physical mecha¬ 
nism of turbulent viscosity can be connected with various 
instabilities. 

The basis for various viscous models of AD is the as¬ 
sumption, that the dynamic viscosity 77 is caused by tur¬ 
bulence of medium and that 77 ~ putlt (rq is the char¬ 
acteristic amplitude of the most large-scale turbulent ve¬ 
locity, i t is the spatial scale, p is the density) (Shakura 
& Sunyaev 1973). In this connection the vital importance 
is given to the research of multimode instabilities, which 
result in a complex perturbations structure. Different spa¬ 
tial and temporal scales may lead to the development of 
turbulence in the disk. 


The thin disk model is the widely used for studies 
of the accretion disk dynamics. The MTD involves av¬ 
eraging over ^-coordinate of the three-dimensional hydro- 
dynamic equations, if a number of additional conditions 
is fulfilled (Shakura & Sunyaev 1973; Gor’kavyi & Frid¬ 
man 1994; Khoperskov & Khrapov 1995). In the con¬ 
text of two-dimensional (in the plane of the disk) mod¬ 
els without magnetic field, four unstable modes of oscil¬ 
lations are present: two acoustic modes (Wallinder 1991; 
Wu & Yang 1994; Khoperskov & Khrapov 1995; Wu et 
al. 1995), a thermal mode and a viscous one (Lightman 
& Eardley 1974; Shakura & Sunyaev 1976; Szuszkiewicz 
1990; Wallinder 1991; Wu & Yang 1994). The instabil¬ 
ity growth rate of acoustic waves increases with the re¬ 
duction of wavelength A. However, the thin disk model 
imposes the restriction on wavelength, A>/i, and there¬ 
fore it is necessary to consider AD 2 -structure for a cor¬ 
rect treatment when A < h. The thin disk model applies 
only to pinching oscillations. Then the perturbed pres¬ 
sure is a symmetric function, and the displacement of gas 
does not move the mass centre in a disk with respect 
to a symmetry plane (2 = 0). Thus, bending oscillations 
(AS-mode) are excluded from consideration, for AS-mode 
the perturbed pressure is antisymmetric function. Also, 
in two-dimensional models the high-frequency (reflective) 
harmonics with characteristic spatial scales in 2 -direction 
< (0.5 4-1 )h cannot be investigated. 

In this paper we investigate the dynamics of acoustic 
perturbations taking into account the non-homogeneous z- 
structure of a viscous disk. Apart from the special problem 
of the limits of MTD applicability, main question is the 
existence of instabilities in short-wave region and, in ad¬ 
dition, the stability of high-frequency harmonics. In Sect. 
2 we define the AD model and we choose the viscosity 
law. In Sect. 3 we consider the dynamics of linear acoustic 
perturbations, and formulate a mathematical problem of 
eigenfrequency determination for various unstable modes 
in the disk. Lastly in Sects. 4 and 5 we discuss results 
of the numerical solution of the boundary problem and 
summarize the main conclusions. 
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2. Model and basic equations 

We shall consider an axisymmetric differentially rotating 
gas disk in the gravitational held of a mass M. Without 
including self-gravity and relativistic effects, and adopting 
cylindrical coordinates we have: 


plane and coincides with a-parameter of the standard the¬ 
ory of the disk accretion (Shakura & Sunyaev 1973). For 
aij it is possible to write down: 
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G is the gravitational constant, ttk = \JGM/r a is Keple- 
rian angular velocity. 

We shall use the axisymmetric hydrodynamic equa¬ 
tions in view of viscosity. The equations of motion and 
continuity have the form 
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where — = —+u—+w—, V = ( u,v,w ) is the velocity, 
dt ot or dz 

P is the pressure, p is the volume density of matter in a 
disk, (jjj is the components of symmetric viscous stress 
tensor (tr ij — <Tji). 

We shall add the thermal equation to the system of 
Eqs. (2)-(5) as 


d S_Q 
dt ~ T ’ 


( 6 ) 


where S is the entropy, T is the temperature, and the 
variable Q defines sources of heat. 


2.1. Equilibrium model and viscosity law 

We shall assume that the equilibrium velocity in the disk 
has only r and <p components: Vq = (Uo, Vo, 0). The ex¬ 
pression for components of the viscous stress tensor may 
be written in the following form: 

Oij = — P , = const > 0 , (7) 


here In = <9(lnO)/<9(lnr). One must be limited by the 
case of weak dependence of equilibrium velocities (Uo and 
Vo) on ^-coordinate for correct transition to the thin disk 
model. This is possible if the following conditions: a rz <C a 
and acp Z <C a are fulfilled. Therefore it is possible not to 
take into account in the equations the terms, containing 
o rz and Ocp Z . The conditions Uq ~ a(h/r) 2 Vo and h/r <C 1 
are carried out for the majority of the models of accretion 
disks. In view of the made above assumptions the equilib¬ 
rium balance of the forces is defined by the equations: 
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The equilibrium functions can be submitted in the 
self-similar form: fo(r,z) = foi{r)fo 2 (z ), if it is assumed 
h/r <C 1. The disk ^-structure depends on the state equa¬ 
tion and the energy flux. For the sake of simplicity we 
shall limit ourselves to a polytropic model 
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Then for the pressure and density it is possible to write 
the solutions: 


Po{r, z) = P 0 {r, 0) F{z) a , p 0 {r, z ) = p 0 (r, 0) F{z) b , (13) 

where F{z) = 1 — z 2 /h 2 , a = n/(n — 1), b = l/(n — 1), n 
is polytropic index, and h defines the disk boundary — in 
the points z = ±/i the equilibrium pressure and density 
are equal to zero. We assume below that h(r) = const. 
The relation 

C 2 (r, z) = ^ = C 2 (r ,0) F(z) , C 2 s (r, 0) = 1 n 2 h 2 , (14) 
Po a 

defines the adiabatic sound speed (7 is adiabatic index). 
For the equilibrium velocities (Vo and Uo) from Eqs. (9), 
(10) taking into account Eqs. (13) and (14) we obtain: 


where i,j = (r, ip, z). As without the account of the second 
(volume) viscosity the trace of the viscous tensor is equal 
to zero (SijOji = 0), thus an = 0 (o rr = o w = a zz = 0). 

The parameters a rz , a vz and a rv> determine a level of 
turbulence in the disk. Moreover a rz and a vz are caused 
by a shear character of flow in ^-direction, and the value 
a r(p is connected to differentiality of rotation in the disk 


Vo (r, z) = rQ(r, z) = r 2 Ct 2 k + IpC 2 ^, (15) 
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3. Linear analysis 


In the framework of standard linear analysis the pressure, 
density and velocity are represented as: 


u = U 0 (r,z) + u(r, z, t ), 


v = V 0 (r, z ) +v(r,z,t ), 
w = w(r, z, t ), 

P = P 0 (r,z) + P(r,z,t ), 

P = Po{r,z) + p(r,z,t). 

In the linear approximation (|/| <C |/o|) from Eqs. (l)-(5) 
the linearized equations become: 
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where g = d^/dz. 

As we shall study the dynamics of acoustic oscillations, 
it is possible to set Q = 0. Then the thermal Eq. (6) 
becomes: 
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here So = cy ln(Po/Po) entropy of equilibrium gas. 
We shall exclude S from Eq. (21) with the help of the 
equation of state S = S(P , p). In the linear approximation 
we obtain: 


S = 
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( 22 ) 


where cy = T{dS/dT ) p , cp = T(dS/dT)p are specific 
heats at constant density and pressure. 

Short-wave approximation in a radial direction ( kr 
1, k is radial wave number) makes it possible to write the 
solution as: 


f(r, z , t) = f(z) expjifcr - i u>t} , (23) 

here to is complex eigenfrequency. 

Taking into account Eqs. (22) and (23), the system 
Eqs. (17)-(21) is reduced to two ordinary differential equa¬ 
tions for amplitudes of the perturbed pressure P(z) and 


the material displacement £( 2 ) from an equilibrium posi¬ 
tion: 
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where Cj = w — kUo, q 2 = to 2 — de 2 , de = fly 7 2(2 + Iq) is the 
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differentiation with respect to ^-coordinate (/' = df /dz), 
£ is the complex amplitude of material ^-displacement 
from equilibrium state, and the following is true w = 
d£/dt = —i u)£. 

Boundary conditions must be added to Eqs. (24) and 
(25). In view of the symmetry, it is natural to consider two 
types of oscillation: 1) the symmetric oscillations (£( 2 ) = 
—£(— 2 ) or P{z) = P(—z)), and therefore 


dP 

«°)=0 or - 


z -0 


= 0, 


(26) 


(such oscillations correspond to the pinch-mode or S- 
rnode); 2) the antisymmetric oscillations (£( 2 ) = £(— 2 ) 
or P(z) = —P(— 2 )), and therefore 


P(0) = 0 or 


d2 


= 0, 
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that correspond to the bending oscillations or AS-mode. 
On the disk unperturbed surface the following condition 
should be fulfilled: 


P(h) + 


z=h 


m = 0 . 


(28) 


By solving the system of Eqs. (24), (25) with the 
boundary conditions Eqs. (28) and (26) (or (27)) (bound¬ 
ary problem), we find the eigenvalues of complex fre¬ 
quency uj for the given distribution of equilibrium param¬ 
eters along the vertical coordinate in the disk. A positive 
imaginary part of the frequency (growth rate) means that 
the eigen-mode is unstable. The above-stated outline will 
have the name ” 3D-model ”, and the model of thin disk 
- ”2D-model ”. 


3.1. Dispersion relation 

For a homogeneous 2 -distribution of the equilibrium quan¬ 
tities, the system of Eqs. (24) and (25) is reduced to the 
following dispersion relation: 

u 4 -La 2 [de 2 +C 2 {k 2 +k 2 z )\-2iujanC 2 s k 2 +de 2 C 2 s k 2 = 0,(29) 
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here k z is the wave number in ^-direction. It is neces¬ 
sary to stress that Eq. (29) at k z = 0 coincides with 
the dispersion relation earlier obtained for 2D-model in 
case of replacement 7 by the flat adiabatic index E = 
r x = (3y — l )/(7 + 1) (Wallinder 1991; Khoperskov & 
Khrapov 1995). Kovalenko and Lukin (1998) have ob¬ 
tained the more accurate formula for the flat adiabatic 
index T = E 2 = (37 - 1 - 7<5)/(7 + 1 — 6) (S = as 2 /n 2 ). 
This relation takes into account Keplerian rotation of a 
thin disk. The adiabatic sound speed in the disk plane 
should be determined by the following expression: 


Ci = r 


Jo P(z)dz 

Jo P( z ) dz 


(30) 


Thus, taking into account Eq. (30), the dispersion rela¬ 
tion Eq. (29) in the case of k z = 0 describes dynamics 
of perturbations within the limits of flat model (MTD). 
Obviously, these oscillations correspond to S-mode. 

In the general case when a > 0, this dispersion rela¬ 
tion gives two unstable acoustic branches of oscillations 
and two damping modes (viscous and thermal ones). The 
damping of the viscous or thermal modes depends on the 
fact that we don’t take into account dissipation and ra¬ 
diative processes in the thermal equation. If these factors 
are taken into account, then the viscous and thermal low- 
frequency oscillatory branches can be unstable (Shakura 
& Sunyaev 1976), but this effect will not change the dis¬ 
persion properties of sound waves. 


4. Results 

It is convenient to characterize the properties of the con¬ 
sidered acoustic oscillations in the disk by the dimension¬ 
less frequency W = ui/fl and the dimensionless wave num¬ 
ber kh. We shall define the basic model as follows: a = 0.2, 
7 = 5/3, h/r = 0.05, l P = -3/2, l p = -1/2, l v = -1/2 
(In = —3/2), n = 5/3. If it is not specified otherwise, the 
parameters take these values. 

4-1. Fundamental S-mode in the 2D- and 3D-models 

Our analysis confirms the good agreement between the 
thin disk model and the results of the 3D-problem solution 
in the case of the dissipational acoustic instability. The 
eigenfrequencies of oscillations obtained from Eq. (29) at 
k z = 0 by taking into account Eq. (30) and from the so¬ 
lution of Eqs. (24), (25) practically coincide in the range 
kh < 3 (see Fig. la,b). The visible differencies occur when 
kh > 3. We emphasize that in this region the formal con¬ 
dition of applicability of the MTD (kh <C 1) is obviously 
broken. In Fig. le,f, the dependencies of group velocity 
v g = R e(du>/dk) and Im (doj/dk) on the wave number k 
for the 3D-model and the 2D-model at T = T x and T = T 2 
are shown. In the case of small values of n, the divergence 


occur at large kh as Fig. la,b shows. The reason is that 
the characteristic scale of inhonrogeneity in the vertical 
direction increases with the reduction of n. The growth 
rate and the phase velocity of perturbations in the frame¬ 
work of the 3D-model is less than for the thin disk model 
(2D-model) in the short-wavelength region. This effect is 
caused by the inhomogeneous distribution of the equilib¬ 
rium quantities in z-direction and by the transverse gravi¬ 
tational force. When the parameter a increases the growth 
rates linearly grow, and the value of the wave number k, 
at which the distinctions between the models appear, does 
not depend on a. The considered low-frequency mode ex¬ 
ists both in the 2D-model and in the 3D-model, and it 
weakly depends on ^-structure. We therefore call it the 
fundamental mode. In the case a«l, the expression for 
frequency is easily obtained from Eq. (29) as: 

, o , ,0^0 . VCfk 2 

u> = ± \Jas 2 + k 2 C^ + ia —~— . 

v s de 2 + C 2 k 2 


This approximation is exact enough. Our results shows, 
that the MTD adequately describes dynamics of pertur¬ 
bations with characteristic spatial scale A ~ 2h 4- r in the 
disk plane. 


The differential rotation (In = — -7 - > 0) and a 

amr 

dependence of dynamic viscosity 77 on thermodynamic pa¬ 
rameters (for example, dependence on density and tem¬ 
perature) are responsible for the instability. A simple 
model can demonstrate this effect. The last term in Eq. 
(3) (i.e., containing a rip ) is a source of the instability of 
acoustic modes. The expression 77 = an (a is the surface 
density, v is the kinematic viscosity) is usually used in 
the 2D-model and for the sake of simplicity we assume 
here that v is constant. Then only a perturbation of sur¬ 
face density a generates perturbation of dynamic viscos¬ 
ity 77 . Now let us consider how a viscous force, which is 
caused by 77 , leads to an amplification of the amplitude 
of sound wave in a disk plane. The evolution equation for 
the surface density, without including Coriolis’s force, is 


d_ 

dt 


d 2 d 2 d 2 a 
dt 2 s dr 2 


d 2 a 

= A——, where A = 2tt 2 lnv 0 > 0. 
dr z 


This equation is linear and the appropriate dispersion re¬ 
lation is u[lo 2 — C 2 k 2 ] = L4/c 2 , which means that acoustic 

1 a 

waves are unstable with the growth rate Im(w) ~ - 7^7 

2 C s 

and entropy oscillations damping (Im(w) ~ —i A/C 2 < 0). 
These estimaties can only be used in the short-wavelength 
limit, because we have neglected epicyclic oscillations. 


4-2. Fundamental AS-mode 

We shall consider low-frequency bending oscillations 
for which Eq. (27) is fulfilled. In a longwavelength limit 
both symmetric and antisymmetric perturbations have 
identical frequency w 2 ~ de 2 = D. 2 in the case of Kep¬ 
lerian disk. This result can easily be obtained from Eq. 
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Fig. la-f. The dependence of the dimensionless eigenfrequency W on the dimensionless wave number kh for different values of 
n (1 — n = 5/3, 2 — n = 1.2) (see a-d). The dependence of dW/hdk on kh for S-mode (see e,f). The solid lines correspond to 
3D-model, dotted lines — 2D-model at T = Ti, long-dashed lines — 2D-model at T = T 2 . 
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Fig. 2a-d. The dependence of W on kh for S- and AS-modes for various values of a. The solid line corresponds to a = 0.2, the 
dotted line — a = 0.1. The numbers of harmonics j are specified in Figure. 


(29) when k = 0. In case of the adiabatic model two other 
branches of oscillations (w 2 = ( C s k z ) 2 ) always damp at 
k > 0 and we do not examine them. The dispersion be¬ 
haviour of S- and AS-modes is very similar (Fig. la-d, solid 
lines). For the fundamental bending mode it is possible to 
set k z h ~ 7t/ 2 in Eq. (29) and for this case the dispersion 
curves on Fig. lc,d are shown as dotted lines. The exact 
solution in the context of 3D-model shows visible differ¬ 
ence from the result obtained from Eq. (29). In the range 
kh ^0.1 these differences are caused by the following fac¬ 
tors: 1) the 3D-model includes radial inhomogeneities of 
equilibrium quantities; 2) the Eq. (29) does not take into 
account vertical inhomogeneity of disk. The obtained re¬ 
sult is very unexpected and significant, since we used Eq. 
(29) beyond the limit of 2D-model approach. 


The physical reason for AS-mode instability is similar 
to the case of S-oscillations. 


4-3. High-frequency S- and AS-modes 

Besides the fundamental S-mode, any number of un¬ 
stable harmonics can be generated in addition. These har¬ 
monics differ from each other by the number of nodes 
of the perturbed pressure across the disk plane. The fol¬ 
lowing estimate of the effective wave number in the z- 
direction can be made: 

k z h ~ 7 rj (S — mode) k z h ~ 7r (j + 1/2) (AS — mode) , 
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Fig. 3a-d. The dependence of W on kh for fundamental mode (j = 0) and for high-frequency one with j = 3 for various values 
of n. Solid line — n = 5/3, dotted line — n = 1.2, dashed line 10. The numbers of harmonics j are specified in Figure. 


where j is the number of the harmonic. The fundamental 
(j = 0 ) and reflective (j > 0 ) harmonics exist for both 
symmetric (S-) and antisymmetric (AS-) modes. 

The dependencies of eigen-frequency u> on the radial 
wave number k for fundamental mode j — 0 and first 
four reflective harmonics j = 1, 2,3,4 are displayed in Fig. 
2a-d. Both the pinch-oscillations and the bending modes 
are unstable (Im(o>) > 0). The imaginary part of the fre¬ 
quency grows with k and reaches a maximum at some 
value. The maximum of Im(w) moves to region of shorter 
wavelengths when the harmonic number j grows, and the 
value increases with reduction of the characteristic scale 
of perturbations in ^-direction. As indicated by Fig. 2a-d, 
the growth rate increases with a, while the perturbations 
phase velocity Re(w)/fc does not depend on a. It should be 


noted that for very short-wave perturbations (kh > 10 ), 
the damping out of oscillations due to the presence of a 
perturbed velocity gradient in the viscous stress tensor can 
be of decisive importance (Khoperskov & Khrapov 1995). 
This factor can be essential for small-scale waves (as on r- 
and on z-coordinate) and it can result in complete stabi¬ 
lization. 

The vertical structure of the disk in the 3D-model 
is determined primarily by parameters 7 and n. In the 
case when n = 7 , entropy does not vary along the , 2 - 
coordinate, i.e. the disk is obviously stable against con¬ 
vective z-motions. When n > 7 we have dSo/dz < 0 and 
the conditions for convective instability are fulfilled, in the 
opposite case n < 7 and dSo/dz > 0, so convection can 
not appear. In Fig. 3a-d, the dispersion curves for various 
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values of n (n > 7 , n < 7 and n = 7 ) for fundamen¬ 
tal modes (j = 0 ) and for the third reflective harmonic 
(j = 3) are shown. The phase velocity of fundamental S- 
and AS-modes weakly depends on the parameter n (see 
Fig. 3b,d), as these perturbations are most longwave in 
^-direction. The adiabatic sound speed in disk C s grows 
with n (see Eq. (14)), and since for the high-frequency re¬ 
flective harmonics k z ~ \di;/i;dz\ yf 0 , the phase velocity 
of perturbations with j > 0 increases also (see Fig. 3b,d). 

The dispersion properties of the acoustic perturbations 
do not depend on the values of parameters Ip and l p in 
the range kh > 0.2. The reason for this is that the radial 
gradients of equilibrium pressure and density give a small 
contribution to the equilibrium balance in the case of a 
thin disk. 


5. Discussion and conclusions 

Acoustic waves in a differentially rotating gaseous disk 
can play an important role for understanding the nature 
of turbulent viscosity in accretion disks. The amplitude 
of small-scale (in r and in z) waves grows most rapidly. 
Such instabilities do not destroy the initial flow at a non¬ 
linear stage, but can effectively make the disk matter tur¬ 
bulent, and in turn, the arising turbulent viscosity gener¬ 
ates unstable sound modes. As a result, a self-consistent 
regime with turbulent viscosity arises. Beside the dissi- 
pational mechanism analysed above, the development of 
global (covering practically the whole disk radially) res¬ 
onant Papaloizou-Pringle modes (Papaloizou & Pringle 
1987; Savonije & Heemskerk 1990), and also resonant am¬ 
plification of acoustic oscillations in the regime of double¬ 
flow accretion (Mustsevoj & the Khoperskov 1991), and 
in the case of disk accretion onto a magnetized compact 
object (Hoperskov et al. 1993) may be important for un¬ 
derstanding turbulent viscosity in the disk. It is remark¬ 
able that in all cases the characteristic time scale does not 
exceed the dynamical time scale (r = 1 /Im(w) ~ H -1 ), 
and that all values are of the same order. 

We have demonstrated the possibility of unstable 
high-frequency acoustic waves in a differentially rotat¬ 
ing gaseous disk. They exist in the system for a limited 
period, and the presence of a positive growth rate does 
not mean that the perturbations must reach a nonlin¬ 
ear stage. The perturbations leave the disk with a speed 
duj/dk ~ C s kh/Vl + k 2 h 2 < C s and the characteristic 
life-time in the disk is ~ r/h ~ 10 2 disk rotation peri¬ 
ods. In view of the derived estimation Im(u;) ~ O.lfl for 
kh > 1 we have the essential growth of wave amplitude 
(exp{0.1r//i|). Nonaxisymmetric perturbations may stay 
in the system for a longer time and reach a non-linear 
stage. 

The following are the main conclusions: 

1. The linear wave dynamics of a thin disk model is com¬ 
pared to an exact solution, which takes into account 


vertical structure. Our analysis shows that not only the 
longwave perturbations, but also rather short waves 
can be considered. A small quantitative difference oc¬ 
cur for a wavelength A < 2irh, and even in the case of 
A > 271r/i the dispersion properties remain qualitatively 
similar. At kh < 1 eigenfrequencies differ less than on 
5%. 

2. In addition to the fundamental dissipational unstable 
sound mode in MTD, an arbitrary number of high- 
frequency unstable harmonics in the case of pinch- 
oscillations is found. These harmonics have a different 
vertical structure. The waves with small scales along 
^-coordinate have the maximum growth rate at short 
wavelengths in the r-direction. The differential rota¬ 
tion and the variable dynamic viscosity cause insta¬ 
bility of all oscillation branches, i.e., a dissipational 
mechanism lead to the perturbation growth. 

3. Taking into account the vertical structure of disk we 
studied new bending modes, which cannot be investi¬ 
gated in the context of MTD. The bending oscillations, 
as well as the pinch-wave, are unstable. The physical 
mechanism of instability and the dispersion properties 
are similar to the case of pinch-oscillations. The signif¬ 
icant feature of all considered unstable acoustic modes 
is the fact that spatial scales of perturbations differ 
from each other, but characteristic growth rates have 
the same order of magnitude. 
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